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ABSTRACT 

Context. More than 10% of extrasolar planets revolve in a binary or multiple stellar system. We investigated the motion of planets 
revolving in binary systems in the frame of the particular case of the three body problem. 

Aims. We analysed of the motion an extrasolar planet revolving in a binary system by following conditions; a) a planet in a binary 
system revolves around one of the components (parent star), b) the distance between the star's components is greater than between 
the parent star and the orbiting planet (ratio of these two distances is a small parameter), c) the mass of the planet is smaller than the 
mass of the star, but is not negligible. 

/Wefhods. The Hamiltonian of the system without short periodic terms was used. Expanded in the terms of the Legendre polynomial 
and truncated after the second order term depending on the one angular variable. In this case the solution of the system was obtained 
and the qualitative analysis of motion was done. We have applied this theory to real extrasolar planets. 

Results. Analysis of the possible regions of motion are presented. It is shown that the case of the stable and unstable motion of the 
extrasolar planets are possible. We applied our calculations to two binary systems hosting an extrasolar planet and calculated the 
possible values for their insufficient orbital elements. For 16 Cyg Bb, there is a region between 44° and 46° for the inclination of the 
planet and the value of the ascending node between 130° and 137°, or the inclination for the planet could have a value from 134° to 
136° and the ascending node from 310° to 317°. For planet HD 19994 b, to have a stable system, the value for the insufficient elements 
should be in the region of 62° to 68° for the inclination of the planet, and the value of the ascending node between 260° and 268°. 
Alternatively, the inclination of the planet also could have the value of 1 12° to 1 18° and the ascending node between 80° and 88°. 
Conclusions. It is shown by which initial Keplerian elements the orbital motion of extrasolar planets could be stable or unstable. 
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1. Introduction 

The motion of an extrasolar planet is calculated using a frame of 
the particular case of the three-body problem. The planet in the 
binary system revolves around one of the components and the 
mass of the planet is much smaller than the mass of the stars, but 
is not negligible. The distance between the star's components is 
greater than between the parent star and the orbiting planet (the 
ratio of these two distances is a small parameter). The motion 
is considered in the Jacobian coordinate system and the invaii- 
able plane is taken as the reference plane. We used the Delaunay 
canonical elements L,, G,, and g, (i=l for the planet's orbit, i=2 
for the star's orbit). They can be expressed through the Keplerian 
elements as 



Lj - Pi yfai , Gi - Li -^1 - e? , Hi - Gi cosli , 
// = Mi , gi = CJi , hi = Q; . (1) 

Where 

mo nil 



Pi ^ k 



y/mo + mi 



(mQ+mi)m2 

= fcyUi , P2= k = kll2 ■ (2) 

yniQ + nil + nil 



In the previous expression the notation has the usual meaning: 
mo, 1112 - the masses of the stars, mi - the mass of the planet, k 
- the Gaussian constant, a,- - the semi-major axis, e,- - the ec- 
centricity, and gi - the argument of the perigee in the invariable 
plane. The eccentricities of the star's and planet's orbits can have 



any value from < e, < 1 . In the general case the motion is de- 
fined by the masses of components and by the six initial values 
of the Keplerian elements of the planet and the distant star. The 
formal solution of a task with the Hamiltonian without short- 
peiiodic terms was obtained in the hyperelliptic integrals by the 
Hamilton-Jacobi method (Orlov, 1988"). Hamiltonian expanded 
in the terms of the Legendre polynomials and truncated after the 
second order terms carries the following form: 

2L\ 2Ll 16 Lpl ^ M ' 

-l5{l-q^){l-T^^)cos(2gi)], (3) 
where the coefficients 71,72, and 73 depend on mass as follows: 



n 



^ ^2 ,2.. ..^2 



Ml 

the parameter: 



72 



Ml 



73 ^k fiiH2-z, (4) 



(5) 



and the cosine of the angle between the plane of the planet orbit 
and the plane of the distant star's orbit 



q = 



2G1G2 



- cos (ii + (2) . 



(6) 
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c is the constant of the angular momentum, and gi is the argu- 
ment of the perigee of the planet orbit in the invariable plane. 
The Hamilton's equation truncated after the second order terms 
depends on the one angular variable only. It is necessary to 
investigate these equations 



dG\ 
dt 



dF 



dgi 
dt 



dF 
'dG[ 



(7) 



In this approximation the task is integrated. This solution we 
have used for intermediate motion. So the semi-major axis ai 
and (22 have no the secular terms; they are restricted. But the 
eccentricity of the planet's orbit can be changed significantly, 
almost as much as 1 . Then the system will be come unstable. The 
change of the eccentricity of an extrasolar planet is described by 
the expression ei - ^\ - ^. Generally, the dependence between 
^ and t is defined by the equation (Orlov J 19881 1: 



-gI 

12 2 



Va 



B3 



1 ym" 
16 



— «! (f- fo). 



(8) 



Where m" = na/ni and ni, M2 are the mean motions of the planet 
and star, ^ = Gj/Lp y = ni2/mo + nii + m2 while, Ai and B3 are 
the constants of integration. 

On the left side of the equation ([8]l under the integral in the de- 
nominator is the square root from the polynomial of the fifth or- 
der. This polynomial can be presented as the product of the two 
polynomials of the second and the third orders, A - fiiOfiiO^ 
where f = \ — e\, and ei is the eccentricity of the extrasolar 
planet's orbit. The establishment of the regions of the motion 
are possible when the roots of the equations /2(^) = 0, f^i^) - 
are found and the signs of the functions defined. 



2. Equations of the second and third order 

We investigated the roots of the equations fii^) - and 
/3(^) = with following forms: 

/2(^) = ^2_2j^2 + 3G^J^ + (c--G2)%^(10+A3)G2 (9) 



and 



Where 



c - 



and 



- Go 
G2 = — 



G2(10+A3) 



€-1(c'-gIJ . (10) 



02 



Pi 



G2 > 1. 



A3 = 2 - 6rilql - 6 (1 - 7/2) [2 - 5 (1 - sir?gio\ 
We noted 



'70 



1 - e 



10 • 



(11) 



(12) 



(13) 



where, eio is the initial value of the eccentricity of the planet, ^0 
is the initial value of the cosine of the mutual inclination between 
the orbits of the planet and it's distant star and giQ is the initial 
value of the argument of the perigee of the planet's orbit in the 



invariable plane. 

So the orbit is assumed to be elliptic, < ei < 1, then < ^ < 1. 
Value of ^ = 1 corresponds on the circular motion. 
For establishment of the boundaries of the regions of the possible 
motion it is necessary to find the roots the equations (|9| and (10 1 



and to define the signs of the function in the interval between the 
roots. 

2.1. The equation of the second order 

We rewrote equation of the second order (j9|l in the following 
form: 



f-2{d^ + 3Gfj^+ -l+2(c2^3G2) 



3Gf\ + \G\h 



.(14) 



Then defined h as difference between values A3 and A3„„, where 



/! = A3 - A 



(15) 



and 

-2 



1 - 8G2 - ^G2q^ - T]^+ 

+20G2{l- q^)sm^gi 



(16) 



We named A^crii the meaning of the constant A3 for which the 
square equation has the root equal to one. This meaning is found 
in the left part of the /2(^) = 0, ^ = 1. We note the roots of the 
equation ( [l4| as €[ and €2, and ei < 62. The coefficient by ^ is 
always less than zero. This equation has no negative roots and 
the free term is also greater then 0. The roots are 



{c^ + 3G2) + ^J{c^ + 3gI + ij 



(17) 



and 62 >> 1 always. Derivation ^ is always positive. So ei is a 

growing function of h. When h - then ei = 1, we conclude for 
h <Q then ei < 1 , and for h > Q then ei > 1 . For the definition of 
the branch of the parabola we found the extreme of the function 

- _2 — 2 - 

in the point which is valid ^ — c + 3G2. The value /2(^) < and 
the branches of this parabola are directed up. 

2.2. Equation of the third order 



We rewrote equation (10 1 using parameter h in the form: 



+ 2P H-G.H- 



/3(^) = f-^{5 + 8c^+4Gfjf 



(18) 



Then multiplied this equation by 4 and identified in its left part 



the item ^/2(^). Then the equation (I81 can be rewrote in the 
following form: 



5(1 -^) 



(c + G2)'-d [^-(c-G2) 



(19) 



For the establishment of the qualitative characteristics of the mo- 
tion, it is necessary to define how the roots of the equation of the 
third order are located in relation to the points of the axes: 



^= 1, 



and 



(20) 
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1. If /z < then we have: 

(a) ^ = 0, /3(0) < 0. 

(b) ^ = eu Hex) > 0. 

(c) f = 1, /3(1) > 0. 

(d) ^ = 62, h{e2) > 0. 

From Sturm theorem (e.g. Dorrie I1965I I it follows that, be- 
tween and €i lay at least one root of the third order equa- 
tion. 

2. If /i > then we have: 

(a) f = , /3(0) < 0. 

(b) ^ = 1 , /3(1) < 0. 

(c) ^ = 61, Me,) < 0. 

(d) ^ - 62 , hie2) > 0. 

The behaviours of the roots of the second and third order equa- 
tions depend on the sign of h and valid: 

1. for /z < is vaHd < 63 < ei < 1 < 64 < 65 < 62. 

2. for /z > is vaHd < 63 < 64 < 1 < ei < 65 < 62. 

We have two roots which are less then 1, therefore we identified 
the roots in ascending order as ^1,^2,^3,^4,^5- This means that 
ei can change in limits from ei„„„ - 1 - ^2 to ei„„„ = 1 - ^i- 
We showed the behaviour of these functions for the hypothetical 
value c, G2 and h on the Figure ( 2.2 1 for h < and on Figure 
(E2I for h>Q. 




Fig. 1. Behaviour of the functions /2(^) (solid hne) and f^i^) 
(dashed line) for h > 0. The left-side graph shows an overview 
while the right-side graph shows the detail between zero and 
two. 
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Fig. 2. Behaviour of the functions /2(^) (solid hne) and fsi^) 
(dashed line) for h < 0. The left-side graph shows an overview 
while the right-side graph shows the detail between zero and 
two. 



3. The investigation of tlie value h 

Consider the value h, which has the following form: 
3(1-7?') 



h 



2G2 

-h20G2(i -^') sin^gi 



1 - 8G2 - 4G2^?7 -1 + 



(21) 



When the value gi 

- _ 3(1 - if) 
2G, 



then value h reaches minimum 



1 - 8G2 - 4G2^?? - 77' 



(22) 



The value of the /!„„„ is always negative. In this case we have the 
two roots valued less then 1, the one root of the equation of the 
third order and one root of the square equation. If - njl then 
h has a maximum value 



3(1-/?') 
2G2 



\-Tf - AGiqil + I2G2 - 2OG2?' 



(23) 



The value of /i„m.v could have both a positive and negative value. 
For the case gi - 7r/2 we rewrote equations (j9]l and ( 10 1 in the 
following form: 



and 

M^) 



(^-'7o)(^-^()-4G2Wo + 8G2) 
+4G2 (1-5^2^(1- 7,2) + I6G 



—2 

2 ' 



(24) 
(25) 



—2 



-3G2^ + G277o?„ (5 - 4^) + 5G2?^ 



(26) 



When one of the roots of the equation (26 1 carries the value of 
^ - tIq, the other roots have to satisfy the equation 

- J j - 770) - 3 G2 ^ + G2 770 qo i5-4O + 5Glql = .(27) 

For hmax the starting value rf - if-^ is one of the boundary limits 
for change in the value of ^. We should establish that ^ - 77^ is 
the least root and the second root has a value less than 1. We 
therefore substitute the value ^ = 77^ in ( 27 1. If the obtained ex- 
pression is negative, then 77^ should become the second root in 



equation ( 26 1, if positive - then it should become the least root. 
The left part of the equation (27 1 for ^ = 77^ is 



5?o - 3770 + ^770^0 (5 - 4?7^) 
G2 



(28) 



Equating this expression to zero and solving the obtained equa- 
tion concerning q^, we find the value of qo of which 77^ is the 



root of expression(28 1. So 



7?0 



q = 



4772 - 5 + 760G2 + (5 - 4772)' 



IOG2 



(29) 



We denote the solution of equation ( 29 1 as qQ\ for the minus sign 
before the root term and the (702 for plus sign. If the value of q^ 
lays inside the region of qm < qo < qoi, expression (28 1 is nega- 
tive, and the value of the root lays between the other two roots. If 
either of the conditions qo < qoi or qo > qo2 is valid, then ^ = 77q 

is the least root of equation (27 1. In this case ei^,^^ = -^1 - 77^ 

and the maximum value of the eccentricity of the planet's orbit 
can not exceed the initial value of eccentricity. The orbit of the 
extrasolar planet may be stable. 
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Table 1. Initial orbital elements of the system 16 Cyg. 

star 16 Cyg A planet 16 Cyg Bb 



Mass [Ms,,,,] 
Mass [Mj„i*smi,] 
Semi-major axis [AU] 
Eccentricity 
Inclination 
Ascending node 
Perigee 
Period 



1.53 

754.53 
0.863 
135.44° 
313.44° 
26.6° 
13512.7 years 



1.68±0.07 
1.68±0.07 
0.689±0.011 



83.4±2.1° 
799.5 days 



Notes. Mass of the parent star 16 Cyg B is 1.01 ± O.O4M5,,,,. We used 
the planet's semi-major axis value 1.6923 AU for the calculations. 



Table 2. Our proposal for a full set of orbital elements for planet 
16 Cyg Bb. 

First possibility Second possibility 



Mass[Mj„p] 

Semi-major axis [AU] 

Eccentricity 

Inclination 

Ascending node 

Perigee 

Period[day] 



2.38±0.04 
1.693 
0.689±0.011 

45°±1° 
133.5°+3.5° 
83.4±2.1° 
799.5 



2.38±0.04 
1.693 
0.689±0.011 

135°±1° 
313.5°±3.5° 
83.4±2.1° 
799.5 



4. Application of the theory on real extrasolar 
planets 

The knowledge of the six pairs of Keplerian elements of the or- 
bits of the system allow us to investigate the character of the 
evolution of the planet's orbit and the possible conditions of sta- 
bility. They may be presented by the orbital parameters, which 
we obtain from the analytical theory. They are - the angle of 
the mutual incUnation between orbits of the planet and star, the 
angular moment of the star and the maximum value of the ec- 
centricity of the planet's orbit. The growth of the eccentricity of 
the orbit could lead to the destruction of the orbit in perigee from 
tidal forces. In the catalogue of extrasolar planets there is a lack 
of data regarding the insufficient elements for the longitude of 
the node and the inclination. Our theory allows for a range of 
possible values for the insufficient elements. In our application 
of the theory, we have selected two binary stellar systems, with a 
planet revolving around one of the components. The first is sys- 
tem 16 Cyg with an interesting planet 16 Cyg Bb, which has a 
perigee close to 90°. The second system is HD 19994. 

4.1. 16 Cyg 

16 Cyg A (HD 186408) and 16 Cyg B (HD 186427) are both 
members of a well-known wide binary system with stars spectral 
types G1.5V and G3V. For the calculations, we used orbital ele- 
ments used by Hauser and Marcy ( 119981 1, which were published 
in the Sixth Catalogue of Orbits of Visual Binary Stars (Hartkopf 
and Mason 2001 ). The semi-major axis was calculated to be 
754.53 AU using the values 21.41 pc (Fuhrmann et. al |I998l l for 
distance and 35.242" for angular separation. We made a revision 
calculation for the semi-major axis using the third Keplerian law 
and the value of a 13512.7 year period (Hartkopf and Mason 
1200 II I. We reached a value of 774.06 AU. There is an obvious 
difficulty in determining orbital parameters for periods greater 
than 10,000 years. Only a small fraction of the orbit has tran- 
spired since astrometric measurements were first obtained in the 
early 1800s. So, we used the value 754.53 AU, derived from the 
distance, for the following calculations because we considered 
the determination period less precise than the determination of 
the values for angular separation and distance. 
The secondary star in the 16 Cyg binary system, 16 Cyg B (par- 
ent star), with a mass of 1.01 ±0.04 Ms,m (Fuhrmann et. al 1998), 
is known to host a giant gas planet 16 Cyg Bb with taxonomy 
class 2J0.2W7 (Plavalova I20I2I I. The extrasolar planet has a 
minimum mass Mjup * sin(i) = 1.68 + 0.07M/„y, in an orbit with 
high eccentricity e = 0.689 (Wittenmyer et. al 2007 ). The start- 
ing conditions for investigating this planet, we took from The 
Extrasolar Planets Encyclopaedia (Schneider et. al 12012b . The 



values of the inclination and ascending note for the planet are 
unknown. In Table ([T]|. are shown the orbital elements for the 
planet and distant star. 

We made a revision calculation for the semi-major axis of the 
planet, using the third Keplerian law and values for the orbital 
period of the planet, which is 799.5 days and published masses 
for the planet and mother star (Table [T]). We calculated it to be 
a I = 1.6923 AU. The mass of the planet is defined as mini- 
mum mass and if we varied this mass to the value of 5 Mjup we 
would get a value of 1.6941 AU. According to the results of pre- 
vious calculations, we decided to use the value 1.6923 AU for 
the semi-major axis and 1.68 Mj„p for the mass of the planet for 
our calculations. 

We varied the elements of the unknown values for 16 Cyg Bb 
in our calculations; the value for the ascending node from 0° to 
360°, 1° at a time, and the value for the inclination from 0° to 
180°, also for each degree. We recognized the pairs of the el- 
ements which when applied to the planet, its motion would be 
unstable, and the others possibly stable. For example for the val- 
ues, 

h =_44°, and Q.i = 134°, _ 
the value h is near to /zmax (h = -0.5696, hi„„x - -0.5698). 
We can use the conditions of stability ^0 < ^01 ilo = -0.9999, 
?oi - -0.5614), and the expression (28 1> and the maxi- 
mum value of the eccentricity is 0.689. With these elements the 
planet's orbit would be stable. Contrariwise, for example, for 
values, 

/i =40°, andQi =20°, _ 
there is a difference between values h and /imax (h = 0.2129, 
hma.x - -0.5721), qo > qoi (qo = -0.3670, qoi - -0.5603), and 
^imm - 0.9295. With these elements the planet's orbit would be 
unstable. 

We found the pairs of the unknown elements, where the one root 
of the second and third order equations have identical values. 
The following are the results we obtained (Table |2]i. We pro- 
posed for a stable orbit, these two possible sets of values for the 
unknown elements; the inclination for the planet could have a 
value between 44° and 46° and the ascending node from 130° 
to 137°, or the inclination of the planet could have a value be- 
tween 134° and 136° and the ascending node from 310° to 317°. 
These values for a stable orbit correspond to the mass of a planet 
equal to 2.38+0.04 Mj„p (for /j = 44° is 2.42 My,,^ and for 
ii = 46° is 2.34 Mj„p) calculated from a catalogue value of 
Mjup * sin(i) - 1.68. We recalculated the planet's semi-major 
axis with this new mass applied through Keplerian law. The new 
value is 1.693 AU. In Table (|2| is shown a proposal for two full 
sets of orbital elements for this planet. 
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Table 3. Initial orbital elements of the system HD 19994. 

star HD 19994 planet HD 19994 b 



Table 4. Our proposal for a full set of orbital elements for planet 
HD 19994 b. 



Mass [Mi„„] 
Mass [My„p*sin ii] 
Semi-major axis[AU] 
Eccentricity 
Inclination 
Ascending node 
Perigee 
Period 



0.37 

151.51 
0.26 
114.1° 
84.13° 
247.74° 
1420 years 



First possibility Second possibility 



1.68 
1.42 
0.3±0.04 



41±8° 
535.7±3.1 days 



Notes. Mass of the parent star HD19994 is 1.34Mj,„,. We used the 
planet's semi-major axis value 1.4273 AU for the calculations. 



4.2. HD 19994 

The binary stellar system HD 19994 (94 Ceti, ADS 2406 AB) 
contains an A component; a yellow-white dwarf with a mass 
of 1.34 Msun and a B component; a red dwarf with a mass of 
0.37 Msun- For the distant star we used Keplerian elements pub- 
lished by Hale ( 1994V To define the semi-major axis we used 
two methods. Firstly, we derived this value using stellar parallax, 
where the distance was 22.38 pc (Schneider l2012l ) and 6.77" for 
the angular separation (Hale 1 19941) and resulted in a semi-major 
axis of 151.51 AU. Secondly, with the application of Keplerian 
law, using the values from Table (|3]l, we calculated the value to 
be 151.37 AU. The difference between these two values is neg- 
ligible we applied the value 151.51 AU in the following calcula- 
tions. 

The planet HD 19994 b was discovered in 2000 (Oueloz ISOOTT l 
and is orbiting an A component. The taxonomy class is 2J0.2G3 
(Plavalova 2012) and the minimum mass is 1.68 Mjup. This 
planet rotates with a semi-major axis of 1.42 AU with quite a 
high eccentricity 0.3 ± 0.04. For our calculations we used the or- 
bital elements published by Mayor et. al. (|20041l. In Table Q the 
initial Kepler orbital elements for the planet and its distant star 
are shown. As with the first system, we revised the calculation 
for the semi-major axis of the planet using Keplerian law. With 
the values listed in Table we valued the planet's semi-major 
axis at a i = 1.4273 AU. If we varied the mass of the planet to 5 
Mjiip we would get a value of 1.4284 AU. We decided to use the 
value 1.4273 AU for the semi-major axis and 1.68 My„p for the 
mass of the planet in our calculations. 

We varied the values for the inclination and ascending node for 
the planet in our calculations, as in the first system. We obtained 
the following results. The stable orbit can have the following 
value of elements; the planet's inclination could have a value 
from 62° to 68° and ascending node between 260° and 268°. 
Alternatively, the inclination for the planet could have a value 
from 112° to 118° and the ascending node between 80° and 88° 
(Table|4|i. These values for a stable orbit correspond to the mass 
of a planet equal to 1.855±0.045 My„p (for ii - 62° is 1.90 Mj^p 
and for /] - 68° is 1.81 Mjup) calculated from a catalogue value 
of Mjup * sin{i) - 1.68. We recalculated the planet's semi-major 
axis with this new mass applied through Keplerian law. The new 
value is 1 .427 AU. In Table Q is shown a proposal for two full 
sets of orbital elements for this planet. 



Mass[Mj„„] 

Semi-major axis [AU] 

Eccentricity 

Inclination 

Ascending node 

Perigee 

Period[day] 



1.855±0.045 
1.427 
0.3±0.04 
65°±3° 
264° +4° 
41°+8 
535.7+3.1 



1.855±0.045 
1.427 
0.3±0.04 

115°±3° 
84° ±4° 
41°+8 
535.7+3.1 



the orbital parameters, which can be calculated from the formu- 
las of the previous sections. There are the angle of the mutual 
inclination between orbits, the angular momentum of the distant 
star and the maximum value of the eccentricity of the planet's 
orbit. When the value of the planet's eccentricity is close to one, 
in a perigee, the tidal forces of the star could destroy the planet. 
We have suggested the possible values for the insufficient ele- 
ments with which the orbit of the planet would remain stable. 
For 16 Cyg Bb there are a two regions of stability. Firstly, from 
44° to 46° for the inclination of the planet and for the ascend- 
ing node, values from 130° to 137°. Secondly, the planetal in- 
clination value could be between 134° and 136° and ascending 
node from 310° to 317°. We proposed the planet's mass to be 
2.38+0.04 Mjup and the value of its semi-major axis equal to 
1.693 AU. For the second system HD19994, the values which 
would suggest a stable motion are in the region of 62° to 68° 
for planetal inclination and 260° to 268° for its ascending node. 
Alternatively, we could suggest 112° to 118° for the inclination 
and 80° to 88° for the ascending node. For the planet's mass 
we proposed the value of 1.855+0.045 Mj^p and for the planet's 
semi-major axis a value of 1.427 AU. 
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5. Conclusion 

We have shown that an extrasolar planet revolving in a binary 
system around one of the components, may rotate on a stable 
or unstable orbits. The conditions for a stable motion depend on 
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